#	  Author: Arturas Rozenas, ar71@duke.edu
#	  Date: Aug 2011
#	  Replication code for paper "Statistical Model for Party 	  System's Analysis", Political Analysis
#    THIS CODE REPLICATES TABLE 2 and FIGURES 2 and 3 in the paper
#    setwd to the current file folder

source("orcomp.R")
source("aux_functions.R")

load("data.japan") # load the data list
					 # contains: 
					 # n - raw number of parties
					 # v - vote-share matrix
					 # X and G - design matrixes for n and y

# the data are from http://www.fps.chuo-u.ac.jp/~sreed/DataPage.html
# they were reshaped and collapsed to suit the analysis: code available upon request
# data were contstrainted to 1960 - 1993 as in Cox 1997
# introducing lags of n, reduced the data set as NA values had to be omitted

library(mnormt)
library(TeachingDemos)

n=as.matrix(data.japan$n)
X=as.matrix(data.japan$X)
v=as.matrix(data.japan$v)
G=as.matrix(data.japan$G)

# NOTE: might take a while to run (more than 7K latent parameters to estimate)

n.sample=2000
ptm <- proc.time()
out=orcomp(n,X,v,G,l.bound=2,n.sample=n.sample,burn=5000,thin=5)
proc.time() - ptm

#========================================================
# run the following to reproduce table 1
#========================================================

D=L=max(n)-1; K=ncol(X); N=nrow(X); l.bound=2

bpost=apply(out$b,2,mean)
gpost=apply(out$g,c(1,2),mean)
mupost=apply(out$mu,2,mean)

cis=function(x){emp.hpd(x,conf=.95)}

conf=signif(round(t(apply(out$b,2,cis)),2))
R1=cbind(signif(round(bpost,2),2),paste("[",conf[,1],", ",conf[,2],"]",sep=""))
rownames(R1)=c(colnames(X),"ln(w)")
R1=R1[c(1,3:5),]
conf=signif(round(t(apply(out$mu,2,cis)),2))
R2=cbind(signif(round(mupost,2),2),paste("[",conf[,1],", ",conf[,2],"]",sep=""))
R2=R2[1:4,]
R=cbind(R1,R2)

n.med=function(b,X){
min(round(median(rtnegbin(10^4,mu=exp(X%*%b[1:K]),dispersion=exp(b[K+1]),l.bound=l.bound))),D+1)
} #### evaluate median predicted count


v.pred=function(b,g,X){ ### predicted vote-share function
n.pred=n.med(b,X)
G.pred=c(X[c(1,3:5)],n.pred)
y.pred=(G.pred%*%g)[1:(n.pred-1)]
return(c(y.v(y.pred),rep(0,L+1-n.pred)))
}

dist=function(x,y){ ##### distance function
return(.5*sum(abs(x-y)))
}

V=matrix(NA,N,L+1) ###### in sample predicted voteshares
for (i in 1:nrow(X)){
V[i,]=v.pred(bpost,gpost,X[i,])
}


print(R)
cat("ln(w)", mean(out$b[,K+1]))
cat("nu",mean(out$nu))
cat("Correctly predicted", mean(abs(V-v)<.05))
cat("MAE", mean(mean(abs(V-v))))

#========================================================
# run the following to reproduce figure 2
#========================================================

m=X[,"m"]
ur=X[,"ur"]
n_t1=X[,"n_t1"]

nlag=function(x,y){median(n_t1[m==x & ur==y])}

X.pred=function(M,UR,n.lag){
c(1,n.lag,M,UR,(UR-mean(ur))*(M-mean(m)),rep(0,10))
}

differ=function(b,g,X.1,X.2){
v.pred(b,g,X.1)-v.pred(b,g,X.2)
}

mn=median(n_t1)
X.3.u=X.pred(3,2,mn)
X.5.u=X.pred(5,2,mn)
X.1.m=X.pred(4,1,mn)
X.4.m=X.pred(4,3,mn)

delta.r=delta.u=matrix(NA,n.sample,L+1)
for (i in 1:n.sample){
delta.r[i,]=-differ(out$b[i,],out$g[,,i],X.3.u,X.5.u)
delta.u[i,]=-differ(out$b[i,],out$g[,,i],X.1.m,X.4.m)
}

M1=t(apply(delta.r,2,cis))
M2=t(apply(delta.u,2,cis))
P=1:(L+1)

m1=apply(delta.r,2,median,na.rm=TRUE)
m2=apply(delta.u,2,median,na.rm=TRUE)


e=.05
plot(seq(0,9,length=L+1),rep(0,L+1),lty=2,type="l",xlim=c(.9,7.3),ylim=c(-.27,.1),xlab="Party rank",ylab="Change in vote-share")
lines(P-e,m1, type="p",pch=4,cex=1.5)
lines(P+e,m2, type="p",pch=1,cex=1.5)
segments(P-e,M1[,1],P-e,M1[,2],lwd=2.2,col="grey60")
segments(P+e,M2[,1],P+e,M2[,2],lwd=2.2)
points(5,-.15, type="p",pch=4,cex=1.5)
text(5.3,-.15,"- M", cex=1.5)
points(5,-.21, type="p",pch=1,cex=1.5)
text(5.3,-.21,"      - Urban", cex=1.5)


#===================================================
#    run the following to reproduce figure 3
#===================================================

v.star=function(m,ur){
t=sample(1:n.sample,1)
return(v.pred(out$b[t,],out$g[,,t],X.pred(m,ur,nlag(m,ur))))
}

e=.1
pd=matrix(NA,n.sample,4)
for (i in 1:nrow(pd)){
pd[i,1]=ifelse(dist(v.star(3,1),v.star(5,1))<e,0,1)
pd[i,2]=ifelse(dist(v.star(3,1),v.star(3,4))<e,0,1)
pd[i,3]=ifelse(dist(v.star(5,1),v.star(5,4))<e,0,1)
pd[i,4]=ifelse(dist(v.star(3,4),v.star(5,4))<e,0,1)
}

C=signif(round(apply(pd,2,mean),2),2)
C

M=c(3,5)
U=c(1,4)
k=1
par(mfrow=c(length(U),length(M)))
par(mar=c(2.5,4.5,2,3))
par(oma=c(.01,.01,.01,.01))
for (i in 1:length(U)){
for(j in 1:length(M)){
vpred=v.pred(bpost,gpost,X.pred(M[j],U[i],nlag(M[j],U[i])))
names(vpred)=c(1:(L+1))
barplot(vpred[1:6], main=ifelse(i==1,paste("M",M[j],sep="="),""),ylim=c(0,.6),ylab=ifelse(i==1 & j==1,"Rural District",ifelse(i==2 & j==1, "Urban District","")),cex.lab=1.5)
}
}

x=1:4
x[1]=expression(paste(pi(epsilon),' = ',"0.08"))
x[2]=expression(paste(pi(epsilon),' = ',"1.00"))
x[3]=expression(paste(pi(epsilon),' = ',"0.49"))
x[4]=expression(paste(pi(epsilon),' = ',"1.00"))

par(fig=c(0,1,0,1), new=TRUE)
plot(1, xlim=c(0,1),ylim=c(0,1), type="n", axes=FALSE, xlab="", ylab="")

text(.4,.9,x[1],cex=1.3)
arrows(.25,.85,.52,.85,lwd=2,code=2,col="grey70")
arrows(.25,.85,.52,.85,lwd=2,code=1,col="grey70")

text(.9,.35,x[2],cex=1.3)
arrows(.82,.2,.82,.5,lwd=2,code=2,col="grey70")
arrows(.82,.2,.82,.5,lwd=2,code=1,col="grey70")

text(.21,.35,x[4],cex=1.3)
arrows(.13,.2,.13,.5,lwd=2,code=2,col="grey70")
arrows(.13,.2,.13,.5,lwd=2,code=1,col="grey70")

text(.4,.17,x[3],cex=1.3)
arrows(.25,.13,.52,.13,lwd=2,code=2,col="grey70")
arrows(.25,.13,.52,.13,lwd=2,code=1,col="grey70")

